using System;
using System.Collections.Generic;

using Atomic.Libraries;
using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Plotting;
using Atomic.Libraries.Plotting.SpacePlotting;
using Atomic.Structures;

namespace Atomic.Phonons
{
	[Serializable]
	public class FrozenPhononDisplacementGenerator
	{
		private FrozenPhononDisplacementGenerator(Structure baseStructure, SpaceVector waveVector, int n1, int n2, int n3, Structure superStructure, double displacement, params Structure[] displacedStructures)
		{
			BaseStructure = baseStructure;
			WaveVector = waveVector;
			N1 = n1;
			N2 = n2;
			N3 = n3;
			SuperStructure = superStructure;
			Displacement = displacement;
			DisplacedStructures = new ReadOnlyList<Structure>(displacedStructures);

			List<Structure> allStructures = new List<Structure>();
			allStructures.Add(SuperStructure);
			allStructures.AddRange(displacedStructures);
			AllStructures = new ReadOnlyList<Structure>(allStructures);
		}

		/// <summary>
		/// Generates superstructure with displaced atoms.
		/// </summary>
		public static FrozenPhononDisplacementGenerator Generate(Structure structure, SpaceVector waveVector)
		{
			return Generate(structure, waveVector, 0.015);
		}

		/// <summary>
		/// Generates superstructure with displaced atoms.
		/// </summary>
		public static FrozenPhononDisplacementGenerator Generate(Structure structure, SpaceVector waveVector, double displacement)
		{
			if (structure.Sites.Count != 1 || structure.Sites[0].Position.PhysicalCoordinates != SpaceVector.Zero)
			{
				throw new NotSupportedException("Only trivial single atom structures are supported currently.");
			}

			SpaceVector a = structure.Lattice.A.PhysicalCoordinates;
			SpaceVector b = structure.Lattice.B.PhysicalCoordinates;
			SpaceVector c = structure.Lattice.C.PhysicalCoordinates;

			int n1, n2, n3;
			Solve(a, waveVector, out n1);
			Solve(b, waveVector, out n2);
			Solve(c, waveVector, out n3);
			Structure s0 = structure.Repeat(n1, n2, n3);

			// Longitudinal and transversal polarizations.
			SpaceVector e1, e2, e3;
			CreatePolarizations(a, b, c, waveVector, displacement, out e1, out e2, out e3);

			// Ready to construct superstructures with displaced atoms.
			Structure s1 = CreateDisplacedStructure(s0, waveVector, e1);
			Structure s2 = CreateDisplacedStructure(s0, waveVector, e2);
			Structure s3 = CreateDisplacedStructure(s0, waveVector, e3);
			Structure s1m = CreateDisplacedStructure(s0, waveVector, -e1);
			Structure s2m = CreateDisplacedStructure(s0, waveVector, -e2);
			Structure s3m = CreateDisplacedStructure(s0, waveVector, -e3);

			return new FrozenPhononDisplacementGenerator(structure, waveVector, n1, n2, n3, s0, displacement, s1, s2, s3, s1m, s2m, s3m);
		}

		private static void Solve(SpaceVector a, SpaceVector q, out int n)
		{
			for (n = 1; n <= 32; n++)
			{
				Complex d = Complex.Exp(Complex.I * SpaceVector.Dot(q, n * a));
				if (Complex.Abs(d - 1.0) < 1.0e-8)
				{
					return;
				}
			}

			throw new Exception("No superstructure solution found.");
		}

		private static void CreatePolarizations(SpaceVector a, SpaceVector b, SpaceVector c, SpaceVector q, double d, out SpaceVector e1, out SpaceVector e2, out SpaceVector e3)
		{
			// Longitudinal polarization.
			e1 = SpaceVector.UnitVector(q) * d;

			// Choose first transversal polarization. Must not be zero.
			e2 = SpaceVector.Cross(e1, a);
			if (SpaceVector.Norm(SpaceVector.Cross(e1, b)) > SpaceVector.Norm(e2))
			{
				e2 = SpaceVector.Cross(e1, b);
			}
			if (SpaceVector.Norm(SpaceVector.Cross(e1, c)) > SpaceVector.Norm(e2))
			{
				e2 = SpaceVector.Cross(e1, c);
			}
			e2 = SpaceVector.UnitVector(e2) * d;

			// Second transversal polarization is easy.
			e3 = SpaceVector.UnitVector(SpaceVector.Cross(e1, e2)) * d;
		}

		private static Structure CreateDisplacedStructure(Structure structure, SpaceVector q, SpaceVector e)
		{
			List<Site> sites = new List<Site>();
			foreach (Site site in structure.Sites)
			{
				SpaceVector delta = e * Complex.Re(Complex.Exp(Complex.I * SpaceVector.Dot(q, site.Position.PhysicalCoordinates)));
				sites.Add(new Site(new DisplacedAtom(site.Atom, delta), new Position(Basis.StandardBasis, site.Position.PhysicalCoordinates + delta)));
			}

			return new Structure(structure.Lattice, sites);
		}
		
		public static IEnumerable<ISpacePlot> PlotLattice(Lattice lattice)
		{
			throw new NotImplementedException();
			/*foreach (SpaceVector p1 in new SpaceVector[] { lattice.A.PhysicalCoordinates, lattice.B.PhysicalCoordinates, lattice.C.PhysicalCoordinates })
			{
				yield return new SpaceLine(SpaceVector.Zero, p1);
				foreach (SpaceVector p2 in new SpaceVector[] { lattice.A.PhysicalCoordinates, lattice.B.PhysicalCoordinates, lattice.C.PhysicalCoordinates })
				{
					if (p1 != p2)
					{
						yield return new SpaceLine(p1, p1 + p2);
						foreach (SpaceVector p3 in new SpaceVector[] { lattice.A.PhysicalCoordinates, lattice.B.PhysicalCoordinates, lattice.C.PhysicalCoordinates })
						{
							if (p1 != p3 && p2 != p3)
							{
								yield return new SpaceLine(p1 + p2, p1 + p2 + p3);
							}
						}
					}
				}
			}*/
		}

		public void Plot()
		{
			/*Structure s0 = BaseStructure;
			Structure ss = SuperStructure;

			// Plot an arbitrary enlarged structure.
			Structure s = s0.Repeat(4, 4, 4);

			// Rebuild longitudinal and transversal polarizations (but larger displacements).
			SpaceVector e1, e2, e3;
			CreatePolarizations(s0.Lattice.A.PhysicalCoordinates, s0.Lattice.B.PhysicalCoordinates, s0.Lattice.C.PhysicalCoordinates, WaveVector, 1.0, out e1, out e2, out e3);

			// Choose one polarization for plotting.
			SpaceVector e = e1;

			List<ISpacePlot> plots = new List<ISpacePlot>();
			if (s0.Lattice is FaceCenteredCubicLattice)
			{
				plots.AddRange(((FaceCenteredCubicLattice)s0.Lattice).PlotWignerSeitzCell());
			}

			// Lattice of base structure and the required superstructure.
			plots.AddRange(PlotLattice(s0.Lattice));
			plots.AddRange(PlotLattice(ss.Lattice));

			SpaceVector lambda = 2.0 * Math.PI * SpaceVector.UnitVector(WaveVector) / SpaceVector.Norm(WaveVector);
			plots.Add(new SpaceArrow(SpaceVector.Zero, lambda));

			// Show atomic positions in a large repeated structure.
			plots.Add(new SpaceDataPlot(new List<Site>(s.Sites).ConvertAll<ISpacePoint>(site => site.Position.PhysicalCoordinates)));

			foreach (Site site in s.Sites)
			{
				SpaceVector delta = e * Complex.Re(Complex.Exp(Complex.I * SpaceVector.Dot(WaveVector, site.Position.PhysicalCoordinates)));
				plots.Add(new SpaceArrow(site.Position.PhysicalCoordinates, site.Position.PhysicalCoordinates + delta));
			}

			Gnuplot.SurfacePlot(plots);*/
		}

		public Structure BaseStructure
		{
			get;
			private set;
		}

		public SpaceVector WaveVector
		{
			get;
			private set;
		}

		public int N1
		{
			get;
			private set;
		}

		public int N2
		{
			get;
			private set;
		}

		public int N3
		{
			get;
			private set;
		}
	
		public Structure SuperStructure
		{
			get;
			private set;
		}
		
		public double Displacement
		{
			get;
			private set;
		}

		public ReadOnlyList<Structure> DisplacedStructures
		{
			get;
			private set;
		}

		public ReadOnlyList<Structure> AllStructures
		{
			get;
			private set;
		}
	}
}
	